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We compute the propagation and scattering of linear gravitational waves off a Schwarzschild 
black hole using a numerical code which solves a generalization of the Zerilli equation to a three 
dimensional cartesian coordinate system. Since the solution to this problem is well understood it 
represents a very good testbed for evaluating our ability to perform three dimensional computations 
of gravitational waves in spacetimes in which a black hole event horizon is present. 
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I. INTRODUCTION 



Biliary black hole systems are among the most promising sources for gravitational wave detectors currently under 
construction such as LIGO, VIRGO and GEO. Theoretical templates of the waveform emitted during the inspiral 

\L} ' and coalescence of such binaries are needed both for increasing the probability of a detection and for extracting 
astrophysical information from the signal |lj . The prediction of such waveforms has therefore become an important 

,__! ■ task of numerical relativity and is the goal of the Binary Black Hole Grand Challenge Alliance 0. 

t^. \ In general, the solution of Einstein's field equations, a large set of coupled nonlinear partial differential equations, 

<3\ i is a task of considerable difficulty. Additional complexity is introduced by the presence of singularities and the black 

hole horizon boundaries necessitated by these singularities. Most significantly perhaps, the computational resource 

limitations, even on large parallel systems, put severe constraints on the accuracy achievable with three-dimensional 

(3D) simulations M. As a result, the problem is far from being solved in general and, at best, there are tailored 

(3JT), solutions to specific difficulties. 

Current successful 3D numerical relativity computations involve either Schwarzschild black holes (i.e. in the absence 
of gravitational radiation) or gravitational waves on spacetimes where no black holes are present. In our calculation we 
combine both components and study the propagation of 3D linearized waves in the fixed background of a Schwarzschild 
black hole. The solution to this problem is well understood; it has been extensively investigated in the past via 
perturbation theory p|-|Io|] and numerous ID calculations have also been made [pll pi. In particular, the gravitational 
waves can be treated as three dimensional perturbations of the background metric and be expanded in terms of tensor 
spherical harmonics. This reduces the problem to solving a ID radial equation for each component of the tensor 
spherical harmonics. For odd parity perturbations this radial equation is known as the Regge- Wheeler equation, 
while for even parity perturbations it is the Zerilli equation. An alternative, non-perturbative approach has employed 
the matching of a Cauchy solution of Einstein's equations onto characteristic hypersurfaces p0[ . 

Here, to test our ability to track gravitational radiation numerically in a 3D black hole spacetime, we (artificially) 
reintroduce the angular dependence into the Zerilli equation and evolve the resulting equation in three dimensions. 
Many of the difficulties which have beset a fully self-consistent, nonlinear calculation are absent in this test problem. 
In particular, the equations are linear and the location of the black hole is known independently of our numerical 
solution. At the same time, our test problem shares many features of the full nonlinear problem. For example, 
we excise the black hole from the computational grid and impose boundary conditions on the apparent horizon. In 
addition, we use the same computational infrastructure usually implemented in present 3D numerical relativity codes, 
which allows both for parallel applications and, in principle, adaptive mesh refinement. Finally, since we can solve 
the ID Zerilli equation with almost arbitrary resolution, we can always compare our 3D results with an essentially 
"analytic" solution. 

The point of this paper is to present this testbed problem |2l]] and demonstrate its potential usefulness as a calibrater 
of numerical accuracy for wave propagation in 3D black hole spacetimes. The numerical exercises we perform are 



illustrative only and our implementation is somewhat crude; also, we do not fully exploit the capabilities of the largest 
parallel machines. Our examples are sufficient, however, to convey the basic idea so that future code builders can 
already begin utilizing this tool in their 3D diagnostic work. 

The paper is organized as follows: in Section [Il| we review the key equations. Section III is devoted to a discussion 
of the different boundary conditions used in the computation, while in Section IV we give a brief description of the 
3D code and of our choices for the numerical implementation. Illustrative numerical tests and results are presented in 
Section M, where we also discuss the relevance of appropriate inner and outer boundary conditions. Finally, Section 
VI contains our conclusions and suggestions for future extensions. 



II. THE ZERILLI EQUATION 

Regge and Wheeler [Q, examining the stability of the Schwarzschild geometry against small nonspherical per- 
turbations, introduced a decomposition of the perturbations in terms of tensor spherical harmonics. Splitting the 
tensor harmonics into "even" and "odd" parity classes, they found that it was possible to write the equations of the 
odd-parity harmonics in terms of a single second order differential equation describing wave motion in an effective 
potential (i.e. the Regge- Wheeler equation). Zerilli || reconsidered the problem and found the corresponding equa- 
tion for even-parity tensor harmonics (i.e. the Zerilli equation). Finally, Moncrief JlO( ] derived versions of both even 
and odd-parity equations using an elegant gauge-invariant formulation of the perturbations. Here we concentrate on 
the even-parity perturbations and describe them in terms of the Zerilli equation, 
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where t and R are the time and radial Schwarzschild coordinates respectively and r* 
coordinate, defined as 



is the so-called "tortoise" 
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The Zerilli function Qi is, in typical applications, constructed from metric perturbations and their derivatives. The 
Zerilli potential Vi(R) is given by 
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where we have used the abbreviations 
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In order to generalize the Zerilli equation to three dimensions, we first replace this tortoise coordinate by the 
Schwarzschild radial coordinate R 
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In order to generalize equation (|g) to a 3D Schwarzschild polar coordinate system (t, R, 9, 4>), we define a new Zerilli 
function Q — Q(t, R, 9, 4>) to be the product of two separable functions 



h{t,R)A{6,<t>) 



(7) 



where Qe(t,R) is a solution of equation M) and A(6,(f>) contains the ("artificial") angular dependence expressed in 
terms of scalar spherical harmonics. We demand that the angular part in the Zerilli function is an eigenfunction of 
the perpendicular differential operator Vj_, i.e. 
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we then write equation m) in terms of 3D Cartesian coordinates 
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In order to take advantage of evolution schemes designed for first-order hyperbolic systems, we rewrite our system 
using the following derived quantities: 



Qo 
0, 



dQ 
dt 
dQ_ 
dx 1 



i= 1,2,3 



(13) 
(14) 



Equation ( |ll|) is then equivalent to the first order system 
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III. BOUNDARY CONDITIONS 



We require three different kinds of boundary conditions. The first ones result from restricting the computational 
domain to one octant, the second ones are the "outer" boundary conditions specified at large distance and the third 
ones are the "inner" boundary conditions on the apparent horizon. We will describe several different choices for the 
inner and outer boundary conditions and compare results in section M. 



A. Octant Symmetry Boundary Conditions 

Octant symmetry permits us to restrict our analysis to a cubical grid in which all coordinate values are positive (this 
is what we refer to as an "octant"). As a result, boundary conditions have to be imposed on the octant-symmetry 
planes and the specific conditions on the functions Qo,Q x ,Qy and Q z can be derived from the symmetry of each 
function on each plane. For even £, Q is symmetric across each plane and hence Qo and derivatives tangential to the 
plane are symmetric as well, while the derivative perpendicular to the plane is antisymmetric. 



B. Outer Boundary Conditions 

At the outer boundaries we impose outgoing wave Sommerfeld conditions. In this approximation, we assume that 
the functions are of the form <f> = fit — R, 9, <j))/R, or equivalently, 

d$ Id 

m + R8R^ = °- < 17 > 

Transforming equation ( J17| ) into a 3D Cartesian coordinate system yields 

m + R x d¥ + R = ° (18) 

and requires, for each grid point on the external faces of the cubical grid, evaluating derivatives perpendicular as 
well as tangential to the surface. However, if / ~ fit — R), then on most parts of the outer surface the tangential 
derivatives are fairly small. It is therefore adequate to neglect these and replace (Eq) with the simpler expression pq] 
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(see Appendix for details on the finite difference form) . 

One concern which might be related with this prescription for the outer boundary conditions is that, strictly 
speaking, it is going to be satisfied only by those grid points on the outer faces which happen to be aligned with 
the radial direction of propagation of the wave-like quantities. All of the other grid points (e.g. those on the edges 
and corners of the cubical grid) will not satisfy expressions (119) identically and a certain amount of reflection will 
occur. While there are several different ways of handling these additional complications [e.g. use of a spherical outer 
boundary embedded within the cubical grid, or an explicit computation of the radial derivative in expression ( |17| ) 
through interpolations f26|| 1, experience has shown that, provided that the outer boundary is placed at a sufficiently 
large distance, the amount of reflection produced is usually very small and conditions ( |l9| ) are sufficient to provide a 
stable outer boundary. 

C. Inner Boundary Conditions 

At the inner boundary we use a horizon excising method which has been discussed by a number of different authors 
[]2T| — 133|| in conjunction with the implementation of apparent horizon boundary conditions. In general, (e.g. as in the 
case of moving black holes) , such a region is not known a priori and its location has to be computed at each time step 
with "apparent horizon finder" routines B4-pq] . This is not the case for the present static Schwarzschild background 
and we need to excise the region inside the event horizon only once during our time integration. Given a masked out 
region of spacetime, the simplest inner boundary conditions involve using suitable finite difference stencils and higher 
order extrapolation methods that fill, on each spacelike hypersurface, the gridpoints just inside the masked region and 
adjacent to the event horizon. These values can then be used in a centered evolution scheme to update the gridpoints 
just outside the horizon. We have implemented such boundary conditions using a fourth order extrapolation scheme. 
These boundary conditions do not violate causality since no information is extracted from within the event horizon 
and, when stable, provide boundary values that are mathematically correct. This prescription is simple to implement, 
does not require special assumptions on the behaviour of the variables in the proximity of the horizon (as would, 
instead, a Sommerfeld condition) and has been proven to be stable for wave propagation in 2+1 dimensions on a flat 
spacetime |n|. 

Alternatively, we can explicitly evaluate equation (O) on the horizon. Choosing Qi = (i = 1,2,3) initially, we 
find the following "freezing" conditions on the horizon 

d t Qo = , 

at R = 2M (20) 

d t Qi = , i=l,2,3. 

These conditions do not involve extrapolation, are very easy to implement into the adopted Macormack evolution 
scheme (see Appendix 0) and as we will show in Section M, produce more reliable results than the more general 
boundary condition based on extrapolations. A different choice of time coordinate may not necessitate such care on 
the horizon. 



IV. NUMERICAL IMPLEMENTATION 

We use different codes to solve the ID and the 3D Zerilli equations. Since the ID code can be run with an essentially 
arbitrary number of gridpoints and hence essentially arbitrary accuracy, we use this solution as an "analytic" solution 
for comparisons. Moreover, since the ID code is considerably simpler, we used this code to experiment with several 
computational techniques before implementing them in the 3D code. 

The 3D code implements a time evolution scheme for equations ( |l5| ) in an environment very similar to other codes 
of the Alliance, in particular the so-called "Empire" code. The latter evolves Einstein's equations in a hyperbolic 
formulation | J22|]23| , so that the mathematical structure of these equations is identical to that of ([H]). Our code 
has been implemented using the DAGH environment p4| , which has been developed for the Alliance. DAGH is a 
Distributed Adaptive Grid Hierarchy software package which allows for parallel applications and, in principle, Adaptive 
Mesh Refinement (AMR). Our code runs in parallel, but we have not implemented AMR here. The code uses a 3D 
cartesian cell-centered uniform grid of typically (32) 3 , (64) 3 , (128) 3 or (256) 3 gridpoints. For most applications we 
restrict the computation to a single octant and a typical run with (128) 3 gridpoints distributed over 16 processors of 
the Origin2000 at NCSA would evolve up to t = 100 M in about 6 hours of CPU time. 

Both the ID and the 3D codes use a Macormack evolution scheme. In this algorithm a "predictor-step" evolves the 
variables, to linear order, from time level t to the subsequent time level t + At and a "corrector-step" uses both the 
old and the predicted values to improve the integration and make it second order accurate (see Appendix |A| for the 
explicit finite difference form of the equations). All of the spatial derivatives during the predictor and the corrector 
steps are one-sided and therefore only first order accurate in space. However, combining the predictor and corrector 
steps on a uniform grid cancels the first order error terms, so that the scheme becomes second order accurate both in 
space and time. 

We verify the second order accuracy of our code in Figure [jL For this particular test we evolved a spherical wave 
on a flat background [i.e. equation ( |l5| ) with c\ — 1 and C2 = c% = 0], for which the solution is known analytically. 
In a second order accurate code the deviation from the analytic solution decreases by a factor of four when the grid 
resolution is doubled. As a test, we plot in Figure u\ the error for three different grid resolutions and multiply the 
errors with successive factors of four. The curves refer to the solution at time t = 3.125 which is reached after 20 time 
steps on the coarsest grid and 160 time steps on the finest. Small deviations between the curves are caused by higher 
order error terms. For increasing resolution, these deviations decrease and the curves converge to a limiting graph. 
This proves that our code is indeed second order accurate. 

V. NUMERICAL RESULTS 

In this section we present numerical results of our 3D evolutions of the Zerilli equation and in particular we compute 
the scattering of gravitational radiation off a Schwarzschild black hole p],||. Consider a gravitational wave packet 
coming from large radius and moving towards the effective potential produced by the black hole. As a result of the 
scattering process, part of the wave packet will be transmitted across the potential and eventually cross the event 
horizon and part of it will be reflected by the potential and propagate out again to spatial infinity. As in quantum 
mechanics, the transmission and reflection coefficients will depend on the wavelength of the incoming radiation (i.e. 
very high frequency modes will be almost completely transmitted while very low frequency modes will be almost 
completely reflected) as well as on properties of the scattering potential (i.e. on the mass of the black hole). The 
knowledge of the rate at which the gravitational radiation is reflected off the black hole as a function of the wavelength 
of the incoming radiation would provide a distant observer with information about the mass of the black hole. 

As initial data for the Zerilli function we choose a Gaussian of width a centered at radius Rq 
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where Pg m is an associated Legendre function. For all calculations in this section we chose Rq = 10M and a = 1M. 
Defining this Gaussian in terms of the tortoise coordinate guarantees that Q vanishes on the horizon. We also assume 
time symmetry, i.e. 

22%£M SOO (0,W,*) = ,22, 

at 

and impose outer boundary conditions at x = y = z = 20 M . 

As the wave packet evolves in time, it splits, with one part of it going towards null infinity and the other reaching 
the horizon, where it induces the quasi-normal ringing of the black hole. The scope of this computation is to calculate 
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FIG. 1. Convergence test of the 3D code. All of the curves show the leading (second) order truncation error for the 
amplitude of a spherical wave Q and of its time derivative Qo- The curves refer to the solution at time t = 3.125 and are the 
result of an interpolation on a fixed set of (16) 3 grid points of solutions computed using (32) 3 , (64) 3 and (128) 3 gridpoints 
respectively. In order to maintain a consistent treatment of the interpolation error introduced, the solutions obtained with 
(16) 3 gridpoints have not been used for this test. 
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FIG. 2. The Zerilli function Q2 as a function of t at radius i? = 15M . The solid curve is the "analytic solution" and the 
others are 3D solutions computed with (32) 3 , (64) 3 and (128) 3 gridpoints. We use one-sided differencing at the inner boundary. 



the waveform of the scattered gravitational radiation as observed at some distance from the black hole and to compare 
it with the "analytic" waveform, i.e. the solution of the ID Zerilli equation ([j]). 

In Figure |] we show Qg with I — 2 as a function of time at radius R = 15 M . For this calculation we used one-sided 
differencing at the inner boundary. More precisely, we use a quartic extrapolation to fill a fictitious gridpoint just 
inside the horizon, so that the gridpoint just inside the horizon has all the neighbors required for a centered updating 
scheme. We can then apply our normal "interior" evolution scheme to update this gridpoint. 

Here and in the following three figures, the solid curve is the "analytic" solution, which we found by integrating 
the ID Zerilli equation (in tortoise coordinates) with a large number of gridpoints. The other curves are the results 
of 3D computations with (32) 3 , (64) 3 and (128) 3 gridpoints. 

The first peak and the first minimum in this curve are produced by the part of the wave packet which moves 
outwards and leaves the computational grid |37|j . The second peak and all the following ones are the quasi- normal 
ringing of the black hole excited by the infalling part of the wave packet. 

Note that for times larger than t « 35M our 3D solutions do not converge to the analytic result. Higher resolution 
improves the solution for early times and delays the time at which the computed curve starts to deviate from the 
analytic one. This behavior is caused by the properties of the Schwarzschild coordinates. Since the "lapse" function 
vanishes on the horizon, infalling waves slow down as they approach the horizon and will, in terms of the Schwarzschild 
time t, never reach the horizon. Stated differently, the characteristic speeds vanish on the horizon, so that the infalling 
waves pile up in front of the horizon. This causes the formation of increasingly small features close to the horizon. 
These will ultimately be smaller than any (uniform) grid resolution. For any constant grid resolution there is therefore 
a time after which features close to the horizon can no longer be resolved. In this region of the spacetime and in 
its domain of dependence, the numerical solution will therefore be spoiled. This is an unavoidable feature of the 
underlying Schwarzschild coordinates together with a uniform grid. 

Obviously, using a high order extrapolation close to the horizon will produce spurious results and will further 
decrease the quality of the solution. In Figure g we show the results of the same calculation, except that we imposed 
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FIG. 3. Same as Figure 2, except that we impose the "freezing" boundary conditions at the inner horizon. 



the "freezing" boundary conditions (|20|), which do not involve extrapolation. As expected, our results improve and 
converge to the analytic solution up to a later time of about t w 45M. 

After this time, the solution can only be improved by dramatically increasing the resolution close to the horizon. In 
these coordinates, the only feasible way to increase the resolution sufficiently, given the memory limitations of today's 
computers, is to use Adaptive Mesh Refinement techniques. We leave to future investigation the application of AMR 
to this problem. 

In Figure we show results, again using the "freezing" boundary conditions, for £ = 4. Since these waves have 
smaller structures and features to begin with, the numerical solution becomes unreliable at an even earlier time. Note, 
however, that these deviations are mostly due to a phase error, while the amplitudes of the waves are reproduced 
quite accurately. 

In addition to a pointwise comparison between analytic and numerical amplitudes, we can compute the energy loss 
through a large sphere inclosing the black hole (Tjjl^]. For £ = 2 this is given by 



dt ~ 384tt V dt 



(23) 



In Figure ra we show the total energy radiated through a sphere of radius R — 15M . Here, our numerical results 
converge to the analytic one even for late times. This is because very little energy is contained in the late-time 
oscillations which have significant phase errors and, later, amplitude errors. 

This becomes apparent in the inset of Figure o, which shows the absolute relative error between the expected 
radiated energy and the computed one as a function of time and grid resolution. Apart from the initial larger values 
related to the outgoing part of the initial wave packet, the relative errors tend to reach constant values of 40%, 10% 
and approximately 2% for (32) 3 , (64) 3 and (128) 3 points respectively. From these estimates and from the fact that 
we are measuring perturbations which have a typical wavelength A ~ 17M, we can infer that a ratio (h/X) ~ f 0~ 2 
on a uniform grid should be sufficient to provide a physical description of gravitational wave propagation on a curved 
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FIG. 4. Same as Figure 3 for £ = 4. 



-4 



-12 




20 40 60 80 100 

t(M) 



-0.4 



0.4 



0.8 

Log Lt(M)J 



1.2 



1.6 



FIG. 5. Logarithm of the energy radiated through a spherical surface at R — 15 M for £ = 2 as a function of time, 
inset shows the absolute relative error between the expected radiated energy and the computed one. 
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background with an error smaller than a few percent. It should be noted that relative errors larger than the ones 
discussed above for the radiated energy can be measured when comparing the computed and "analytic" waveforms at 
a given time and location. However, these errors appear only at later times, when the amplitude of the perturbation 
is drastically reduced and its contribution to the total radiated energy is negligible. 

Besides allowing us a somewhat more accurate simulation, for any given resolution, the freezing boundary conditions 
( P0| ) have also solved the instability problems encountered with the use of one-sided differencing at the inner boundary. 
With "freezing" boundary conditions we were able to evolve the code up to a physical time t ~ 2 x 10 M which 
corresponds to rs 10 5 crossing times. At that stage the L^ norms of the relevant variables were showing a negative 
slope, clear indication of stability of the code. 

A final comment should be made about some modifications to the present approach which could improve the quality 
of the numerical simulations. As discussed in Sections O and M, one of the major difficulties encountered in this study 
are related to the choice of a static slicing of the Schwarzschild background. A first improvement can come from 
adopting a coordinate gauge that better describes wave propagation in the vicinity of the horizon and that avoids the 
"freezing" of the constant coordinate-time slices (see, e.g., p8|). With a harmonic background, the wave perturbations 
will propagate through the horizon in finite coordinate time rather than piling up there. Although the wave equation 
becomes more complicated, it should be numerically straightforward to repeat our test problem in this context. A 
second independent improvement could come from a more ingenious use of the memory resources. We have shown 
that the use of high resolution grids improves the agreement between the computed and the analytic solution. AMR 
techniques are well known to be particularly suitable to study those configurations in which a very high resolution is 
necessary only in some regions of the computational domain. It is likely that their implementation in the proximity of 
the horizon would provide the resolution necessary to adequately resolve many of the details of the wave propagation. 
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VI. CONCLUSION 

We have presented 3D computations of gravitational radiation scattering off a Schwarzschild black hole. This test 
problem represents a useful tool to investigate numerical issues such as finite differencing across the horizon, inner 
and outer boundary conditions, evolution schemes, code parallelization, optimal use of numerical resources and could 
be used as a standard benchmark on 3D numerical relativity codes. The numerical results obtained from this code 
are in good agreement with the "analytical" ones and converge to the latter as the resolution is increased. Late time 
deviations from the analytic solution are due to specifics of our choice of coordinates. 

The present study has also allowed us to test the minimum computational requirements of 3D numerical relativity 
computations against the present resources available on modern parallel computers. We have found that simple 
physical configurations such as a perturbed Schwarzschild black hole can be handled quite satisfactorily with most 
of the physical content of the solution being reproduced. However, we also believe that much greater resources than 
the ones available today may be necessary in order to study more complex perturbative problems and, of course, the 
fully relativistic evolution of binary black holes. 
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APPENDIX A: FINITE DIFFERENCE EQUATIONS 

We here give explicit expressions for the finite difference forms implemented in the second order accurate Macormack 
evolution scheme. We use a standard notation where lower indices refer to the spatial location of the gridpoint and the 
upper index refers to the time level. As mentioned in Section fill we solve for a set of 5 first order partial differential 
equations [i.e. (JTa)] that have the symbolic form 

dB _dA dC _0A dDdA 

dt ~ dx ~ ° ' dt ~ dy ~ ° ' at ~ dz ~ ° ■ (A2) 

All variables are first evolved from the time level t = t n to the time level t = i n+ i = t n + At by means of the predictor 
step, in which the spatial derivatives are computed using a first order accurate backward differencing: 



A Ul = A ti,k + At [^ B ^ ~ S * n -u,J + ^(Cta - Ch-xj.) + ^- z (D?,j, k - D^ k _ x ) + (A3) 



KRS(A i: j k ,B i: j k , Cij^f.) > , 

Bi,j,k = Bi,j,k + ~^~ \^H,j,h ~ ^i-l,j,k) j (A4) 

/nrra+l _ (~m , f±^_ I in _ 411 \ ( \K\ 

D t,j,k = D?,j,k + ~^1 {^i,j,k ~ A?,j,k-l) ■ ( A6 ) 

New, second order accurate values of the variables at the time level t = t n+ i are then computed with the corrector 
step, in which first order forward differencing is used for the spatial derivatives: 
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The first order error term introduced in the predictor step is exactly cancelled in the corrector step, so that the 
algorithm becomes second order acc urat e. Th e fif th equation of ( |l5| ) involves a simple time integration; its finite 
difference form can be deduced from (A3) and (A7) when all of the spatial derivatives are taken to be zero. 



APPENDIX B: OUTER BOUNDARY CONDITIONS 



We here present second order accurate finite difference expressions for the outgoing-wave Sommerfeld conditions of 
the form (M) . In general, these conditions need to be applied at the six external faces of the cubical grid, but here 
we will concentrate only on the expressions for the (x Max , y, z) plane, from which equivalent expressions for the other 
outer planes can be derived. 

After some manipulations involving a correct centering in time and space of the relevant variables, the outer 
boundary conditions (H9) for the gridpoints (i M , j, k) on the (x Max , y, z) plane assume the form 
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where 



and 



(R) 
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: ,k -Ri M -l,j,k 



F= (^)At 
(a:) h 



(x) 



C i M J. fe 



-l,j,k 



(B2) 



(B3) 



In the case in which "octant" symmetries are not used, the equivalent expressions for the gridpoints (i m ,j,k) on the 
(x mtn ,y, z) plane will be obtained by changing (i M , j, k) -> (i m , j, k), (i M - 1, j, k) -> (i m + l,j, k) and F -► -F. (i M 
and i m are here the first and last grid points in the ^-direction respectively.) 
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